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Abstract 



This work aims to provide standard formulations for direct minimization approaches on various types of static problems of contin- 
uum mechanics. Particularly, form-finding problems of tension structures are discussed in the first half and the large deformation 
problems of continuum bodies are discussed in the last half. In the first half, as the standards of iterative direct minimization 
strategies, two types of simple recursive methods are presented, namely the two-term method and the three-term method. The dual 
estimate is also introduced as a powerful means of involving equally constraint conditions into minimization problems. As exam- 
ples of direct minimization approaches on usual engineering issues, some form finding problems of tension structures which can be 
solved by the presented strategies are illustrated. Additionally, it is pointed out that while the two-term method sometimes becomes 
useless, the three-term method always provides remarkable rate of global convergence efficiency. Then, to show the potential ability 
of the three-term method, in the last part of this work, some principle of virtual works which usually appear in the continuum me- 
chanics are approximated and discretized in a common manner, which are suitable to be solved by the three-term method. Finally, 
— — l some large deformation analyses of continuum bodies which can be solved by the three-term method are presented. 
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1. Introduction 

Within this work, standard formulations for solving vari- 
ous types of static problems of continuum bodies by the direct 
minimization methods are presented. The direct minimization 
methods are always associated with static mechanics via prin- 
ciple of virtual work. For example, the direct minimization ap- 
proaches are sometimes very effective on solving form-finding 
problems of tension structures] 1]. In particular, the aim of this 
work is to present the basic strategies such as the three-term 
method and the dual estimate, and to illustrate various types 
of static problems that can be solved by using them. 

In section 2, as the standard recursive direct minimization 
methods, the two-term method and the three-term method 
are described. While the former is basically identical with the 
steepest decent method and the latter is with the dynamic re- 
laxation method, some differences are pointed out. In addition, 
via discussion of a form-finding problem of a simple cable-net 
structure as a typical example, the relation over the principle 
of virtual work, the stationary condition, and the standard 
search direction is clarified. Furthermore, the dual estimate 
is proposed as a powerful means of involving constraint condi- 
tions into direct minimization approaches. Then, form-finding 
analyses of a tensegrity structure and a tensioned membrane 
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structure are illustrated as examples of minimization problems 
with constraint conditions. 

In section 3, more general cases of static problems of con- 
tinuum bodies are taken into account. First, the discrete prin- 
ciple of virtual work, the stationary condition, and the stan- 
dard search direction are formulated as the result of standard 
procedures. They can be positioned as the generalizations of 
those appeared in section 2 and enables the direct minimization 
methods feasible on general cases of static problems of contin- 
uum bodies. Finally, some large deformation analyses of con- 
tinuum bodies which can be solved by the three-term method 
are illustrated. 

2. Two-term method and three-term method 

2.1. Direct minimization approaches without constraint condi- 
tions 

Suppose a form-finding problem of a prestressed cable-net 
structure which can be stabilized via introducing prestress (see 
Fig. 12. U . For example, any solutions of the following station- 
ary problem of a functional can be used as such a form: 



n(*) = £ 



WjL) (x) 



stationary, 



(2.1) 



where wj, Lj denote the weight coefficient and the length of the 
j-th cable respectively. The weight coefficients are the param- 
eters which are assigned with the aim of varying the form by 
varying them, and they are treated as constant in the following 
formulations. In addition, x is a column vector which contains 
the unknown variables {x\, ■ ■ ■ ,x„}. 
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In this work, x and corresponding gradient vector are al- 
ways arranged as 



x = [ X\ ■ ■ ■ x„ J , and V/ = [ 



ILL 

d.\\ 



' £]■ (2-2) 

By the authors, it has been pointed out [ 1 ] that solving Eg. ( 12. Il l 
by the direct minimization methods can be positioned as that an 
equilibrium equation provided by the force density method Q 
is represented in a different manner firstly, and then the equilib- 
rium equation is solved by a direct minimization method which 
differ from the method proposed in the original force density 
method. 

In this work, the unknown variables \x\ , • • • , x„] are always 
assumed as denoting the Cartesian coordinates of the free nodes. 
In addition, remark that those of the fixed nodes are eliminated 
beforehand from x and they are directly substituted into each 

When Eq.( I2. lb is solved by the direct minimization meth- 
ods, n (x) is usually called the objective function. Additionally, 
the direction of greatest rate of increase of II, namely 



_ Vn 

r ~ W\ 

is usually adopted as the standard search direction. 

The stationary condition of Eq. ti2.lb is as follows: 



(2.3) 



vn = o o Yj 2w i L F L i = 



(2.4) 



Here, taking the inner product of Eq. ( 12.41 ) with arbitrary col- 
umn vector Sx, namely 

Sx = [ 6x\ ■■■ 6x„ \ , (2.5) 
the principle of virtual work can be obtained as: 

5w = 2wjLj6Lj = 0, (2.6) 



or the variational principle can be found has 



where 



sn = o, 



Sf = Vf ■ Sx 



(2.7) 



(2.8) 



is the variation of /. Due to the arbitrariness of Sx, Eq. (|27< 
and Eq. ( 12.7b are always equivalent with Eq. 12.41 It is amazing 
that the common frameworks which are provided by the clas- 
sic mechanics, such as the principle of virtual work and the 
variational principle, can be even found in such a minor force 
density method. 

On the other hand, the principle of virtual work for self- 
equilibrium cable-net structures can be expressed as 



Sw = TijSLj = 0, 



(2.9) 



where nj denotes the tension of y'-th cable. By comparing Eq. 
( 12.6b and Eq. ( 12.9b . when II (x) is stationary, at least one self- 
equilibrium state can be found as 



n m J = [ 2w\L\ ■■■ 2w m L m J. (2.10) 



Therefore, any solution of Eq. ( 12.1b can be used as a form of 
cable-net structures that can be prestressed. 

When the standard search direction is given by Eq. (12.3b . 
one of the simplest recursive direct minimization methods is 
given by 



vn 7 



° (X — JCCurrent) i 



ivni 

-^Next — -^Current — Q^Current) 



(2.11) 



which is called the two-term method in this work. Here, "Cur- 
rent" and "Next" are the current and the next step numbers. As 
is immediately noticed, the two-term method is basically iden- 
tical with the steepest decent method. The main differences are 
as follows: 

• The standard search direction is always normalized. 

• Step-size factor a is a parameter which is assigned with 
the aims of adjusting the rate of convergence and treated 
as constant in the formulations. (In the steepest decent 
method, step-size is usually determined by line-search al- 
gorithm ) 

The aim of the normalization of the standard search direc- 
tion is to prevent the divergence of the computation. Moreover, 
without using some computational algorithms to determine a , 
if a is treated as constant in the formulation and to be adjusted 
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by somebody via GUI, appropriate a can be found easily. Ac- 
tually, it was really easy and intuitive operation to determine a 
via GUI. 

By the way, the rate of global convergence efficiency of the 
two-term method is not basically good, as of the steepest de- 
cent method usually is. Because it is supposed that the compu- 
tation would usually start from a point which places far from 
the exact solution, the rate of global convergence efficiency 
must be improved. Then, the following remedy of the two- 
term method sometimes provides a remarkable improvement 
of global convergence efficiency: 



VIT 



'"Current — | VII| ° ~~ * t ~ UIrent ^ ' 
?Next = °- 98 ?Current ~ ^Current, 
-^Next — X Current + a 9Next' 



(2.12) 



which is called the three-term method in this work. When 
\x,q,r] are thought as {position, velocity, acceleration}, Eq. 
( 12.121 ) can be positioned as one kind of equation of motion 
with a damping term, therefore the basic idea of the three- 
term method is almost identical with the dynamic relaxation 
methodJ3l. However, as same as in the two-term method, 
the standard search direction is also normalized in the three- 
term method. Then, it is better to interpret the three-term 
method as just one of the recursive direct minimization meth- 
ods and is not being based on dynamic mechanics. The factor 
0.98 which can be found on the second line means that 2% of 
q is compulsory cut in each step, which can be interpreted as 
one kind of damping factor. This factor is having no basis and 
being determined by some experience. 

As is mentioned above, because the three-term method is 
not based on the dynamics, the following consideration is not 
precise; however, the high rate of global convergence efficiency 
provided by the three-term method can be understood intu- 
itively when it is explained with terms of energy conservation 
law. Namely, due to the elimination of 2% of q in each step, the 
total energy of the system is compulsory exhausted gradually, 
then II and q would shortly leach to the minimum value and 0. 

An numerical model for verification of the two-term and 
the three-term method is provided by Fig. 12.11 which consists 
of 5 fixed nodes and 220 tension members. The coordinates of 
the fixed nodes are also presented in the figure. As shown by 
Fig- 12. 1T b). initial values of {x\,- ■■ ,x n ] were set by random 
numbers ranging from -2.5 to 2.5 by the authors, then it was 
able to obtain Fig. 12. li e) and (d)by either 2-term or 3-term 
method. Fig. 12.11 (c) is the form taking minimum value of the 
sum of squared length of all the tension members and the cor- 
responding minimal value was 160.214. Fig. 12.11 (d) is the form 
which was obtained when 4 times greater weight coefficients 
were assigned onto the boundary cables and the corresponding 
minimum value was 188.09. In both method, 0.2 was used as 
the step-size factor a. 

Fig. 12.21 shows the history of the objective function when 
Fig. 12.11 (c) was obtained. As shown by Fig. 12.21 after a while 
a was fixed to 0.2, soon II converged and vibrated around 160. 



At that time, the norm of |VII| was 0.13, which my be thought 
as not being sufficiently small. Even such cases, as shown in 
Fig. 12.31 it is possible to decrease |VI1| gradually by decreas- 
ing a gradually. However, this work expects the two-term and 
three-term method to be used as a means of exploring various 
equilibrium forms by varying the parameters such as weight co- 
efficients and the coordinates of the fixed nodes freely, in which 
a would be kept constant such as 0.2. 



Figure 2.2: History of U(a = 0.2) (by 3-term method) 




Figure 2.3: History of | VH| (by 3-term method) 



2.2. 



History of 3-term method 

Here, it must be noted the close relation between Eq. ( 12.121) 
and the "Three-term recursion formulae". In 1959, "Three term 
recursion formulaeAh was firstly presented by M. Engeli, H. 
Rutishauser et. al [3]. In 1982, M. Papadrakakis stated that 
the dynamic relaxation method |7[] and the conjugate gradient 
method, can be classified under the family methods with three- 
term recursion formulae [4]. Because Eq. ( 12.121 ) has a com- 
mon form with the conjugate gradient method and its basic idea 
highly resemble the one of the dynamic relaxation method, it 
may be possible to position the three-term method proposed 
in this work as the simplest method based on the three-term re- 
cursion formulae 

2.3. Direct minimization approaches with constraint conditions 





(a) Connectinos (b) Result 

Figure 2.4: Form-finding of Simplex Tensegrity 
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(b) Orthogonal decomposition of search direction 

Figure 2.5 : Direct Minimization Approaches under Constraint Conditions 
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Figure 2.6: Form-Finding of Complicate Tensegrities 



In this section, direct minimization approaches with con- 
straint conditions are discussed. As an example, let us con- 
sider the form-finding problem of a Simplex Tensegrity struc- 
ture, which is shown by Fig. 12.41 A Simplex Tensegrity is a 
self-equilibrium structure that consist of 3 compression mem- 
bers which are shown as thick lines in the figure and 9 tension 
members, as thin lines in the figure. In addition, remark that the 
members are pin-jointed on only their ends. 

In general, it is expected to obtain such a self-equilibrium 
form when some objective function with respect to the lengths 
of the tension members is minimized by constraining the lengths 
of compression members. Then, let us consider the follow- 
ing simple minimization problem with equally constraint con- 
ditions: 



II w (x) 



s.t. 



/ t WjLj 4 (x) — > min, 

7=1 

(Lio-Iio) =0, 
(Ln-L n ) = 0, 



(2.13) 



where \L\, ■ ■ ■ , L9} denote the lengths of the tension members 
and {Lio, • • • , £12} denote the lengths of the compression mem- 
bers. In addition, {w\,- •■ , W9} are the weight coefficients as- 
signed to every tension members. Moreover, \L\q, • • • , Z42} are 
the constraint values of the lengths of the compression mem- 
bers, which are treated as constant in the formulations below 
but are assigned with the aims of varying the form by varying 
them. The basis of the power 4 which is put on Lj is not ex- 
plained in this work, because it has been already reported by 
the authors yj . 
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Fig. 12.41 (b) shows the form that is taking the minimum 
value of Yjj L 4 j when every lengths of the compression members 
are constrained to 10.0. The corresponding minimum value was 
18000 and the methods performed by the authors are described 
below. 

Applying the Lagrange multiplier method, Eq. ( 12.131 ) re- 
duces to the following stationary problem of a functional: 



m+k 

(x) - Lm+k'j -» stationary, 

;=1 k=l 

(2.14) 

where the first sum is taken for all the tension members and the 
second sum is taken for all the compression members. In ad- 
dition, A is a row vector containing the multipliers [A\ , • • • , A r }. 
From now on, let x, A and the corresponding gradient vectors 
be arranged as 



x = [ x\ 
X = \ h 



x„ 



,1, 



ox ~ I aii 1 



. (2.16) 



Then, let the gradient operator V be defined by V/ = |£. 

The stationary condition of Eq. ( 12.14b can be expressed as 
a set of two conditions: 



^ = & ^ = 0. 
ox oA 



(2.17) 



Thus, in this work, it is important that unknown variables x and 
multipliers A are explicitly distinguished and the corresponding 
stationary conditions are discussed separately. 

First, let us discuss the first stationary condition, namely 

Vn = 0«J| 4wjLj 3 VLj + Yj A kVL k = 0, (2.18) 

j k 

and its general form is expressed as 

vn = vn w + a ■ j x = o, (2.19) 

where VII W is the gradient of the objective function and J \\ is 
an Jacobian matrix given by 



J A = 



(2.20) 



which must be refreshed in each step. 

In this work, the number of the constraint conditions are as- 
sumed as being smaller than the number of the unknown vari- 
ables. In addition, bad conditioned problems, such that satis- 
faction of the constraint conditions is almost impossible, are 
not discussed. The simple interpretation of above statements in 
terms of mathematics is that J j is always supposed as full-rank 
and the number of the columns is supposed as being greater 
than the one of the rows. 

Taking into account the making use of the direct minimiza- 
tion methods on this problem, the indeterminacy of VII must 
be solved, i.e. due to the unknown multipliers [Aw - , A r ] , 



Vn can not be determined uniquely, and this makes the two- 
term and the three-term method infeasible. In contrast, if 
an additional rule is supplemented with the aim of determin- 
ing [Ai, ■ ■ ■ , A r ] uniquely, VII is also determined uniquely and 
both the two-term and the three-term method turn to feasible. 
One of the simplest ideas is making use of the Moore-Penrose 
type pseudo inverse matrix J\. 

First, Eq. (12.19t is transformed into 



a ■ j A = -vn H 

and then, A can be determined by 

a = -vn„, • j\ 



(2.21) 



(2.22) 



which provides basically a least norm solution. When J ' \ is sup- 
posed as fullrank, it is simply given by J\ = j\ ■ (j a ■ /JQ . 
One may feel it is very hard to adopt such a least squared solu- 
tion because it is not an exact solution; however, when x turns 
to a solution, A given by Eq. ( 12.22t turns to a least norm solu- 
tion and when Eq. ( 12.22b gives a least norm solution, it implies 
that a stationary point has been obtained. Otherwise, when Eq. 
(12.221 1 gives a least squared solution, it implies that x still has 
not leached to a stationary point, therefore, a supplement of an 
additional rule to determine A uniquely must not be interfered 
by any reason. 

As the result of above discussion, a unique mapping from x 
to Vn can be defined by 



vn = vn w + a-j a = 0o(a = -vn w ■ /+) , 



(2.23) 



which determines a gradient vector filed and thus both the two- 
term and the three-term method turn to feasible. The deter- 
mination of {Vn, A} by using Eq. ( 12.231 1 is essentially identical 
with the dual estimate, which is defined in linear programming 
theory, particularly in the context of the primal affine scaling 
method |0]. 

By the way, the substitution appeared in Eq. (12.23b can be 
performed immediately, and then Eq. (12.231 1 reduces to 



vn = vn w (i-Ja 



Ja + ), 



(2.24) 



which is widely known as the projected gradient in terms of 
the projected gradient method and in which A is eliminated. 
However, the multipliers are always calculated explicitly in this 
work because the dual estimate can be interpreted to the com- 
position of forces when VII,,. is considered as a force, A as a 
reaction force and VII as a resultant force as shown in Fig. 12.51 
(a). 

Let us recall and discuss the second stationary condition, 
namely 

(im+l — ^m+lj = 

: , (2.25) 

\Lm+r ~ L m +rj — 

which is apparently the prescribed equally constraint conditions 
themselves. One of the simplest ideas to satisfy Eq. (12.251 1 is to 
solve simultaneous linear equations such as 



5 



J i ■ Ax = -r, (2.26) 

where Ax is a correction vector of x and r is a residual vector 
given by 

Lm+l (x) — L m+ i 

: . (2.27) 

The definition of J j is apparently identical with Eq. (12.20b . 
but should be refreshed again. Here, the Moore-Penrose type 
pseudo inverse ]\ plays an important role again to determine 
Ax as 



Ax 



(2.28) 



which basically gives a least norm solution. In addition, be- 
cause x can places far from the hyper-surface on which the 
constraint conditions are satisfied, it is highly recommended to 
rescale Ajc to prevent the computation being unstable, such as 



XCu 



:= x Cu 



+ 0.5 Ax, 



(2.29) 



where : symbol represents a substitution of the right hand side 
into left hand side. If Eq. ( 12.291 ) is always performed once af- 
ter the execution of Eq. ( 12.111 ) or Eq. ( 12.12b in each step, x 
would gradually approaches to the hyper-surface on which the 
prescribed constraint conditions are satisfied and soon, the mo- 
tion of x generated by the two-term or the three-term method 
will be constrained onto such a hyper-surface, as shown by Fig. 
I2.5f b). By using either two-term or three-term method, Fig. 
12.41 (b) was obtained. By introducing Eq. (12.291 ), it is also en- 
abled starting the computation from random numbers. As same 
as in the previous section, to obtain Fig. 12.41 (b), the authors 
gave random numbers ranging from -2.5 to 2.5 to the initial 
values of {xi, • • • , x„) and set the step-size factor a as 0.2. 

By comparing Eq. (12.241 ) and Eq. (12.281 ). one may notice 
that Vfl and Ax are row and column vectors which are selected 
from completely decomposed two spaces that are orthogonal to 
each other, because \I — J\ ■ J a) represents the kernel of J\ 
and vice versa. As is depicted in Fig. 12.51 (b). the 71-dimensional 
search space (usually assumed as an Euclidean space) that x 
belongs to is firstly decomposed into a group of hyper-surfaces 
on which residual vector r taking the same value. Then, on 
each point of each hyper-surface, the vector space attached to 
each point is completely decomposed into the tangent subspace 
and the orthogonal complement. Finally, each of Vfl and Ax 
is selected from each of the tangent subspace and the orthogo- 
nal complement respectively. Thus, the feature of the method 
proposed in this work which should be emphasized is that two 
subspaces that are orthogonal to each other correspond to two 
separated stationary conditions, and, two different strategies are 
performed independently on each subspace. 

By the way, when minimization problems with constraint 
conditions are solved by the method proposed above, the two- 
term method become sometimes useless, particularly on com- 
plicate problems. In contrast, by using the three-term method, 
it was still possible to find the forms of complicate structures 
such as the tensegrities shown by Fig. I2.6l (see U]). 



Such complicate problem on which the three-term method 
works better than the two-term method can be easily found 
widely. Fig. 12.71 shows such another form-finding analysis, 
in which, the analytical model consists of cables, membranes, 
compression members and fixed points, and based on the fa- 
mous Tanzbrunnen in Cologne AiFrei Otto , 1959Aj- The se- 
lected stationary problem that was solved is as follows: 



n (x, A) = £ Wj L 4 j (x) + J] w k S 2 k (x) 

7=1 * 
r 

+ X A ' ( Lm+I ^ ~ ^ m+ ') ~* stationar y> ( 2 - 3 °) 

i=\ 

where the cables, the membranes are subdivided into line ele- 
ments and triangle elements, and the first sum is taken for all 
the line elements and the second sum is taken for all the tri- 
angle elements. In addition, the third sum is taken for all the 
compression members, which can be composed into a station- 
ary problem by applying the Lagrange multiplier method to the 
constraint conditions. Moreover, Lj, S k , and L m +i are the func- 
tions that respectively represent the length of a line element, 
the area of a triangle element, and the length of a compression 
member. 

The statinoary condition with respect to x can be expressed 

as : 



vn = J] 4u ; /.;v/. ; + J] 2w k s k vs k + J] a,vl„ 1+i = o, 

j k I 

(2.31) 

and then, taking the inner-product between 6x and Eq. ( 12.31I ). 
the principle of virtual work for this problem can be expressed 
as: 

m r 

5w = Yj ^ w i L ) 5L i + Yu 2w * S * 6S * + Y W*u = 0. (2.32) 

j k 1 

As means of solving Eq. ( 12.30b . while two-term method was 
completely useless, the three-term method worked really fine. 
What is more important is that it was also possible to vary the 
form by varying the weight coefficients or the lengths of the 
compression members and to explore the possible self-equilibrium 
forms. 

On the basis of above considerations, this research is strongly 
focused on the three-term method, even though the two-term 
method, or the steepest decent method, is sometimes described 
as one of the most standard direct minimization methods. In the 
next section, the formulations for solving various types of static 
problems of continuum bodies by the three-term method are 
presented. 



3. Continuum mechanics 

When the gradient of the volume of a tetrahedron element, 
Wj, is added to a set of VLy and VS j, a compact framework 
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in which a set of {VLy, VSj, VV,j is adopted as basic gradient 
vectors can be formed; however, this framework is almost use- 
less, because, while the stress tensor defined on a 3-dimensional 
body usually has 6 degree of freedom, the degree of freedom 
that Wj has is at most only the scalar multiplication. Hence, 
a further consideration on the continuum mechanics must be 
needed, if taking into account the three-term method on solv- 
ing various types of static problems of continuum mechanics. 
In this section, the discrete principle of virtual work, the sta- 
tionary condition, the standard search direction are derived 
from the principle of virtual work that is a governing field 
equation. Additionally, a general form of {VL,, VSj, VV)J is 

presented as oj N j \T]A. 

3.1. Minimal surfaces and uniform stress surfaces 

From now on, Einstein summation convention is used. In 
this subsection, the relation between minimal surfaces and the 
uniform stress surfaces are discussed. A minimal surface is a 
surface such that when its form is varied arbitrarily by fixing its 
boundary, its surface area does not change. 

In general, the surface area of a surface having a fixed bound- 
ary can be expressed as: 

a = ^"da, da = ^det g u d6 l d6 2 (1 < i,j < 2), (3.1) 

where gjj, 8 l , 9 2 represents the Riemannian metric and the local 
coordinate parameters which are defined on each point of the 
surface. Then, the variation of surface area can be expressed 
as: 

Sa = | sJdetgijdO^dd 2 . (3.2) 

J a 

Here, it is widely known that S -y/detgy can be expanded as 



S^detgij = -g u 6g u , 



where g'i is defined as the inverse of namely g'j - (gijj . 
Additionally, Sgij is not completely arbitrary but must be geo- 
metrically admissible, namely dgij must be expressed as 

S gij = (ViSu k ) g kj + [Vj5u k ) g ki - 2h u Su\ (3.4) 

where \Su l , 5u 2 , <5m 3 } are arbitrary scalar fields such that each 
5u k satisfies 5u k = on the boundary; however the detail of 
Eq. (13.41 ) is not discussed in this work because what is only 
needed for the three-term method is just an approximation of 
Sgij and not Sgij itself. While Eq. ( 13.4b is a field equation, the 
aim of this work is to avoid such complicate and difficult field 
equations and to obtain an approximated solution easily by the 
direct minimization methods. 

Substituting Eq. ( 13.3b into Eq. ( 13.21 i. the minimal surface 
problem can be expressed as: 



O 



U/ Ss « 



da = 0. 



(3.5) 



Eq. (13.5b has a close relation with the principle of virtual work 

of self-equilibrium membranes whose boundary is fixed: 



Sw 



da = (1 < i,j < 2), 



(3.6) 



where t and cr' J respectively denote the thickness and the Cauchy 
stress tensor defined on each point of the surface. Additionally, 
\Sg,j is used instead of the variation of strain due to the essen- 
tial identity between them. 

Using raising and lowering indices law of tensors, i.e. X'i = 
X' k g k j, the principle of virtual work is transformed into: 



Sw 



\fy. k g kj 6g if 



da = (1 <i,j,k<2). (3.7) 



(3.3) 



Moreover, when a new stress tensor T' k is defined by 



n = tcr[ k , 



(3.8) 
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Eq. ( 13.71 ) is transformed again into: 



9) 



6w z 



5w' 



On the other hand, Eq. ( 13.5b . the minimal surface problem is 

also transformed into: 



6a = i J^/^da = (1 < j, j, Jt < 2) , (3. 



10) 



which can be a simple demonstration of the essential identity 
between minimal surfaces and uniform stress surfaces. 

3.2. Principle of virtual work for N -Dimensional Riemannian 
manifolds 

In this subsection, the formulations appeared in the previous 
subsection are generalized into any dimensional spaces from 
2-dimensional spaces (surfaces). The length, the area, and the 
volume of a curve, a surface, and a body which have a boundary 
are expressed as: 



l - JtJ^Sgt/ii = (1 < i, j,k < 2) , (3.19) 

\ fotf'Sgijfo = (1 < i,j,k < 3), (3.20) 

where t and A respectively denote the sectional area of a cable 
and the thickness of a membrane. 

Here, when new stress tensor T\ is defined for each dimen- 
sion individually as: 

T\ = Aa-[ k (N = 1) , T\ = to-\ k (N = 2), and T\ = a[ k (N = 3) , 

(3.21) 

Eq. (ED , Eq. (T3~T7b and Eq. ( l3T20b are unified into: 



c5w n 



4 f ^Sgtjto" = (1 < i, j, k < AO , (3.22) 



which is the principle of virtual work for self-equilibrium N- 
dimensional Riemannian manifold M. 

Here, Eq. (13. 17b . the minimal volume problem of M, can 

be transformed into: 



/= fdl, a= Ida, vh f dv, (3.11) Sv" = - f S i k g iJ Sg ij dv N = (1 < i, j, k < N) . (3.23) 

J I Ja J v 2 



where dl, da, dv are respectively called the line element, the sur- 
face element, and the volume element which are defined by 

dl = VHT^ 1 , (3.12) 

da = yjdetgijde^e 1 (1 < i, j<2), (3.13) 

dv = y jdetg u d0 1 d8 2 d0 i (1 < i,j < 3), (3.14) 

where g (/ represent the Riemannian metrices defined on each 
point of each geometry. Such geometries on which Rieman- 
nian metric is defined on each point, can be classified as 1,2,3- 
dimensional Riemannian manifold. 

Noticing the common forms appeared in Eq. (13. lit . (13 . 1 2b . 
( 13.13b . (13.14b , it is very natural to define the volume element 
and the volumen of a A^-dimensional Riemannian manifold M 
by 

dv^ = Jdetgijdd 1 ...d6P, v N s f dv^, (3.15) 
' Jm 



M 

By comparing Eq. ( 13.22b and Eq. ( 13.231 ). it can be noticed that 
the minimal volume problem is a special cases of the princi- 
ple of virtual work such that T\ = 6' k , and the principle of 
virtual work is one of the natural generalizations of the min- 
imal volume problem. In general, 8', can be classed with the 
unit matrix. 

3.3. Galerkin method 

The principle of virtual work which is defined in the pre- 
vious sub section, i.e. 



5w • 



\ f Vte«dv* = o (3.24) 

1 Jm 



then the variation of the volume of M can be expressed by 



5v N = - 



idv". 



(3.16) 



M 



Hence, the minimal volume problem of M can be defined by: 



= \ f S i] Sgif 
L Jm 



6^ = - gv'dgijdy" = (1 < i, j < AO . (3.17) 

i M 



is basically a field equation; namely the degree of freedom of 
8gij is infinite. Then, with the aim of solving the principle 
of virtual work by the direct minimization methods, in this 
subsection, discrete principle of virtual work is deduced. 

First, when the form is explicitly represented by n indepen- 
dent parameters such as [x\, ■ ■ ■ ,x n ), then x,6x, andV/ = -J- 
can be defined with the same manner of section 2. When, the 
degree of freedom of the form is «, then at most n independent 
dgij can satisfy Eq. ( 13.24b . Thus, in general, any form on which 
n independent dgij can satisfy Eq. ( 13.241 ) is usually adopted as 
an approximated solution. One of such natural ways of giving 
6gjj is altering dgij into 



Sgij = Vgij ■ Sx, 



(3.25) 



By the way, the self-equilibrium equations of cables, mem- 
branes, and 3-dimensional bodies whose boundary is fixed can 
be expressed in the form of the principle of virtual work as 
follows: 



which is essentially the Galerkin method. If 6gij is altered into 
Eq. ( 13.251 ), the discrete principle of virtual work (weak form) 
is obtained as: 



6w 



. _ 1 



l - jA^Sgijdl = (i,j,k= 1), (3.18) 6wN = \f T -k8 kj {^8u ■ 6x ) d v N = 0, 



(3.26) 



then, letting 6x out of the integral operator, discrete principle 
of virtual work (strong form) is obtained as: 



[jj M Ti k S k ^8ijdv N] j-dx = 0, (3.27) 



finally, due to the arbitrariness of 6x, 



L Jm 



' k g kj V 8i jdv" = 0, 



(3.28) 



which is the discrete stationary condition and can be also po- 
sitioned as a discrete form of a self-equilibrium equation. 

When external forces are acting on the manifold M, the 
discrete principle of virtual work (strong form) is firstly ex- 
pressed as 



and it follows 



5 L T!k8kiV8ijdyN ) ' 5x = p5x ' 

L JM 



O to 



(3.29) 



(3.30) 



where p is a row vector containing the components of the nodal 
loads, which should be basically derived via some discretization 
process of continuum load but further detail is not discussed in 
this work because it has been already discussed in the usual 
finite element formulations. 

Then, since the discrete stationary condition is an n-order 
simultaneous non-linear equations and the number of the un- 
known variables is n so that basically it can be solved. In ad- 
dition, when the discrete stationary condition is solved by the 
direct minimization methods, 



CO 



(3.31) 



is adopted as the standard search direction. 



3.4. N -dimensional Simplex elements 

The discrete stationary condition to = which is derived in 
the previous subsection still contains integral operator, which is 
the last obstacle to be overcome. In this subsection, as a power- 
ful means of calculating to on general numerical environment, 
Af-dimensional Simplex element is presented. 

When the integral domain is subdivided into m elements, if 
element integral is defined by 

"1 s \ f T a r g^V ga/s dv N , (3.32) 

where the integral operation is calculated separately within each 
element, then, to can be simply expressed as 



to ■ 



The most simplest idea to calculate Eq. ( 13.32b is to let the inte- 
grated function constant within each element. Fig. ( 13. U shows 



1,2,3-dimensional Simplex elements, which are apparently just 
the line, triangle, and tetrahedron element having N + 1 nodes. 

From now on, first, let \p lt ■ ■ ■ ,p N+ \] be a set of the Carte- 
sian coordinates of the nodes of an element, and then second, 
let {tf 1 , • ■ ■ , be a simple local coordinate defined on the ele- 
ment. Third, let each coordinate parameter be taking the value 
from to 1 . Then, finally, the global coordinate (assumed as 
the Cartesian coordinate) of each point within the element can 
be given by an interpolation function defined by 

r(e 1 ,--- ,e N ) = e 1 ( Pl -p 2 ) + --- + e N (p N -p N+l ) + p N+l . (3.34) 
Then, referring to the definition of the base vectors, namely 



= 

g l - ■■ g N can be calculated by 

gi=Pi-P M (l<i<N), 



(3.35) 



(3.36) 



which is apparently constant within the element. Hence, the 
Riemannian metric 

gij = gi ■ gj (3.37) 

is also constant within the element. Moreover, when consider- 
ing the usual elastic bodies, T', is usually dependent on only 
gij, then T\ is also constant within the element. As the result of 
above considerations, the integrated function is constant within 
the element and the following formulations can be used: 

<o)(T a ^ l -Lj[T\g^gn}., (3.38) 

M ) ( r r) s \ S i [TW^g a p]j (1 < a,p, r < 2) , (3.39) 

oj) (r;) = I Vj [T^Vgap]. (1 < a,fi, 7 < 3) , (3.40) 

where Lj,S j, Vj respectively denote the length, the area, and the 
volume of each dimensional element, namely they are given by 



Vdet^ 



Vj = i A /deT^J 



(1 < a,p< 2), 



(1 < a,0<3). 



(3.41) 
(3.42) 
(3.43) 



In addition, g'-i is the inverse of gij. The inverses of tiny matri- 
ces can be calculated by using the following explicit represen- 
tations: 



(3.44) 
(3.45) 







1 












gu ' 


511 512 


1 _ 1 


g22 ~g\2 


521 522 


detgy 


~g2\ g\\ 
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(a)7V=l (b)N=2 (c) N=3 

Figure 3.1: Simplex Elements 



-1 



£11 


gl2 


gli 


1 




gl2 




£13 


£21 
£31 


822 
g32 


g23 
g33 . 


detgy 




g22 
g32 . 


X 


£23 
£33 . 



£13 




£11 




£11 




£12 




£23 


X 


£21 




£21 


X 


£22 




. £33 . 




£31 




£31 




£32 





(3.47) 

where 



a\ 




bi 




aib-i 


- b 2 a 3 




a 2 


X 


b 2 




aj,b\ 


- £>3fli 


(3.48) 


ai 




b 3 




a\bi 


-b\Q2 





Thus, the tiny inverses have been completely eliminated from 

tO. 

3.5. Gradient vectors and the general form 




(a) to '(T r 5<) (b) cofiV^d') (c) o> ( T <) ) 
Figure 3.2: {r k = # k ) 

In this subsection, the relation between to N . (X") and the 
gradient vectors are discussed. 

Interestingly, as are depicted in Fig. 13.21 when T" = 6", the 
following exact relations are formed: 

to) (s° y ) = VLj, (3.49) 

a? } fa) = VS Jt (3.50) 

o](6 a r ) = VVj. (3.51) 

The demonstrations of above relations can be obtained by alter- 
ing 6 symbols into V symbols in the demonstration of that the 
minimal volume problem is a special case of the principle of 
virtual work such that T" = 6 a y , which was described in the 

subsection 3.2. Therefore, a set of IVLj, VSj, VV)} coincides 
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with {co) (T«) , coj (T a y ) , ai* (t° 7 )} when T" = 6° y , and then 

to N j [T"y\ is one of the natural generalizations of {VL ; , VS j, VV ; }. 

Furthermore, to N j \T?y) can be used when {VLy, V5 j, VV,j are 
calculated. 

Based on above considerations, it can be noticed that only 
special cases such that T" is given as just a scalar multiple of 
6 a y have been discussed in the section 2. Therefore, it is very 
natural to consider general functions as T a y . Particularly, a map 
from gjj to T" is no other than the constitutive law it self. 

The explicit representations of VLj and VS j are presented 
in Appendix A and one may notice that they look very different 

while the difference between to 1 , and to 2 is just the dimensions 

j j J 

of the matrices; however, due to g'i which is defined as inverse 
matrix, the apparent difference between VLj and VS j is resulted 
from the difference between the explicit representations of tiny 
inverse matrices. 

By the way, each cay is just a mixture of the gradient vec- 
tors Vgij, then even though a function fj such that a/! = Vfj, 
is not found in usual, to N . is a row vector that highly resemble 
the gradient vectors. Therefore, the discrete stationary condi- 
tion which was formulated in the subsection 3.3 is expected to 
be solved by the two-term or the three-term method by just 
altering VII into to. 

3.6. Numerical examples 

In this subsection, some examples that discrete stationary 
conditions can be solved by the three-term method are illus- 
trated. As the simplest constitutive law, only 

T' k = EgHgik-gik), (3.52) 
:.r u -r,g ki = Eg u (gik-g,k)g k} , (3.53) 

is considered, which is the uniform linear material with Poisson 
ratio=0, and where E is the stiffness factor. It must be remarked 
that the stiffness factor E is identical with Young's modulus 
only when N = 3, otherwise it is multiplied with the sectional 
area or the thickness, when a Riemannian manifold is related 
with a real material. Additionally, e& - {gn - gn) is no other 
than the strain tensor itself and gn is the Riemannian metric 
treated as constant and is measured on the initial shape on which 



the stress tensor vanishes. Note that while T' k is not a symmetric 
matrix, T' is a symmetric matrix. 

Unlike the numerical examples described in the section 2, in 
each initial step of the following numerical examples, {xi , • • • , x„} 
were not given by random numbers but were given by the co- 
ordinates of the initial shape, and g% were calculated on such 
initial shapes. 

Fig. 13.31 and [J4l show natural forms of handkerchief that 
are hanged by 1, or 2 point under gravity. The dimension of 
the numerical model is 8.0x8.0, and every z-components of the 
nodal forces were set as 0. 1 . Each form has been obtained by 
solving 

co = 2^(n)-p = 0, (3.54) 

i 

which is a discrete stationary condition or a discrete form of 
equilibrium equation, by the three-term method. In addition, 
E — 50 as the stiffness factor and a = 0.2 as the step-size factor 
were used. 




(a) Initial Shape (b) Result 



Figure 3.3: Natural Forms of Handkerchief 1 




(a) Initial Shape (b) Result 1 (c) Result2 



Figure 3.4: Natural Forms of Handkerchief 2 

Fig. [33] shows the large deformations of a cantilever under 
gravity whose dimension is 2.0-2.0-12.0. Fig. 13.61 shows the 
large deformations after buckling of a bar which has the same 
dimension of the former . Each form has been obtained by solv- 
ing 

a> = Y J <o 3 j (r k )- P = o, (3 - 55) 

which is the discrete stationary condition, by the three-term 

method. In both analyses, E — 50 as the stiffness factor and 
a = 0.2 as the step-size factor were used. 



In the analysis which resulted Fig. 13.51 every z-components 
of the nodal forces were set as p, which is shown in the fig- 
ure. In the analysis which resulted Fig. 13.61 small random num- 
bers were firstly supplemented to the initial nodal coordinates to 
make the model easily buckle. Then, z-components of the nodal 
forces of only 9 nodes which place on the top of the model were 
set as p, which is shown in the figure. Even if this can be ex- 
plained as one kind of buckling phenomena, the analysis itself 
is just a large deformation analysis; hence precise identify of 
critical load is almost impossible. However, the Euler buckling 
load corresponding to this example was calculated as p cr =1.14 
and its division by 9 is 0.126, which indeed places between Fig. 
[16](a) and (b). 




(a)p=0.0 (b)p=0.05 (c)p=0.1 



Figure 3.5: Large deformations of a cantilever under gravity 




(a)p=0.10 (b)p=0.20 (c)p=1.0 

Figure 3.6: Large deformation after buckling 



4. Conclusions 

In the first half of this work, the direct minimization ap- 
proaches were discussed, in which some form finding problems 
of tension structures were considered. Especially, as the stan- 
dard strategies of for direct minimization approaches, the two- 
term method, the three-term method, and the dual estimate 
were presented. In addition, the relation over the principle 
of virtual work, the stationary condition, and the standard 
search direction were clarified, which are the means of direct 
minimization approaches. 

In the last half of this work, starting from the principle 
of virtual work (field equation) that usually appeared in the 
continuum mechanics, the discrete principle of virtual work 
was deduced. Moreover, the the discrete stationary condition 
and the standard search direction were formulated to let the 
three-term method feasible. Those formulae were expressed 
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with a/j \T'\ which is one of the natural generalizations of 

{VL,, VSj, VVjj, hence, the last half of this work was a gener- 
alization of the first half of this work. Finally, some large de- 
formation analyses of continuum bodies were illustrated. Those 
various types of numerical examples which were shown in this 
work imply the potential ability of the three-term method that 
can be a powerful means of solving various types of static prob- 
lems of continuum bodies. 



and its visualization is presented by Fig. IA.1I 
Let us investigate dL, i.e. 



5L = VL- 



Sp 
6q 



(A.6) 



As shown in Fig. IA.2I 5p and dq are firstly projected to the line 
determined by p and q, then, 6L is measured on the line. 
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Appendix A. Gradients 

Appendix A.l. Gradient of Linear Element Length 
Suppose p and q denote two nodes. Let 





Px 






p = 


Py 


and q = 


q y 




. Pz . 




. 1z . 



(A.l) 



represent the Cartesian coordinates of p and q. 

The length of the line determined by p and q is given by 



L {Px, Py, Pz, <7.v, q y , q : ) 
= ^(Px - q x ) 2 + (p y - q y ) 2 + (Pz ~ <7z) 2 - 
If the gradient of L is defined by 

VL 



1 dL 


dL 


dL 0L dL dL 


'[dp. 


Spy 


dp z ' dq x ' dq y ' dq z 



its components are as follows: 

VL = [ — - — — - % Pz - - — - — — - — q ' - — 1 



(A.2) 
(A.3) 

(A.4) 
(A.5) 




Figure A.l: VLj 




Figure A.2: Variation of Element Length 

Appendix A.2. Gradient of Triangular Element Area 
Let p, q, and r be three vertices. Let 





Px ' 




q x ' 




r x 


p = 


Py 


. q = 


q y 


, r = 


r y 




. Pz _ 




. 9z _ 




r z 



(A.l) 



denote the Cartesian coordinates of p, q, and r. 

The area of the triangle determined by p, q, and r is given 

by 



S(p x ,--- ,r z )=^^N~N, 

(N ={q-p)x{r- p)) . 

If the gradient of S is defined by 



V7C - T i£_ SS_ dS_ 

Vd = [ dp x < dp y > dp z 



ds_ 

dr- 



(A.8) 
(A.9) 



(A. 10) 
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its components are as follows: 



which can be expanded as 



vs = 





1 














(r- g) X I 







1 







1 


































,(p-r)xj 







1 







) 














1 




















,(q-p)X 







1 







1 




o 









1 





where n is defined by 



n = 



N_ 



and a visualization of VS is presented by Fig. IA.3I 
Let us investigate SS , i.e. 



SS = -n-((r- 



q) X 6p + (p - r) X 6q + (q - p) X Sr) . 



(A. 11) 



(A. 12) 



(A. 13) 



With respect to 5p, for example, when dp is orthogonal to the 
element, (r - q)x6p becomes orthogonal to «, then 6S vanishes 
(see Fig. IA.4b . On the other hand, when dp is parallel to the op- 
posite side, (r - q)x dp vanishes, then 6S vanishes. Therefore, 
only the component of 5p which is parallel to the perpendicular 
line from p to the opposite side can produce 6S . In other words, 
5S is measured on the plane determined by p, q, and r. 



dsn = 

dx M (x J+ , - Xj) - dXi (x J+1 - Xj) + dXj+, qc m - x,y - dXj (x i+ , - x ; ) 

+ dY M (y ;+1 - Tj) - dY, (Y 1+l - Yj) + dY j+i (X M - Y t ) - dYj <J M - Y,) 

+ dZ M (Z jtl -Zj)-dZi(Zj^ -Z,)+rfZ J+1 (Z, + , -Z,)-dZj(Z M -Z,), (A.15) 

where {Xi, Yj,Zj} (1 < i < N + 1) represents the Cartesian coor- 
dinates of z'-th node. 

When the independent parameters [x\, ■ ■ ■ , x n ) are selected 
as the Cartesian coordinates of all the free nodes {Xj, Yj, Z,} 
(1 < i < N + 1), by comparing Eq. (IA. 15b and the following 
relation, the explicit representation of Vgn can be obtained. 



d Sij 



d 8U 
dx n 







dx\ 




= Vgij ■ 











(A.16) 




Figure A.3: VS; 



r-q 




Figure A. 4: Variation of Element Area 



Appendix A.3. Gradient of Riemannian Metrics 

The explicit representation of Wg/j can be obtained by refer- 
ring the following calculations. First, Eq. ( 13.351 ) and Eq. ( 13.361 ) 
follows 

dgij = d (p M - p t ) ■ (p J+l - pj) + (p M - pt) ■ d(p J+l - Pj ) , (1 <i,j<N) (A. 14) 
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